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Abstract 


This  thesis  documents  a  new  approach  to  investigate  the  gamma  radiation  activity 
of  the  fission  products  of  three  different  fuels  (U-235,  U-238,  and  U-239)  exposed  to 
three  different  incident  neutron  energy  spectra  (thermal,  fast  spectrum,  and  high 
energies).  An  application  of  the  exponential  moments  function  is  used  with  a 
transmutation  matrix  in  the  calculation  of  complex  radioactive  decay  chains  to  achieve 
greater  precision  than  can  be  attained  through  current  methods.  The  result  of  this 
research  is  a  code  which  can  calculate  the  decay  products  from  complex  radioactive 
decay  chains  with  a  high  degree  of  precision  while  quantifying  the  uncertainty  in  gamma 
activity  due  to  uncertainties  in  the  isotope  properties. 
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PRECISE  CALCULATION  OF  COMPLEX  RADIOACTIVE  DECAY  CHAINS 


I.  Introduction 

The  activity  of  radioactive  isotopes  is  based  on  two  factors:  the  quantity  of  the 
isotope  present  and  the  half  life  of  the  isotope.  This  thesis  attempts  to  document  a 
methodology  for  improving  the  calculated  activity  through  increased  precision  of  the 
calculation  of  each  isotope  at  any  given  time  and  direct  calculation  of  the  radioactivity. 
Also,  the  propagation  of  uncertainty  in  the  activity  due  to  the  uncertainties  in  the  isotope 
data  is  quantified  through  a  Monte  Carlo  method. 

The  half  lives  of  isotopes  span  many  orders  of  magnitude,  from  short-lived 
isotopes  that  decay  nearly  instantaneously  to  stable  isotopes  with  very  long  half  lives. 
Because  of  the  large  range  of  half  lives,  calculating  the  quantity  of  each  isotope  present  in 
a  given  decay  chain  after  a  period  of  time  may  be  difficult  due  to  the  stiffness  the 
differential  equations  describing  radioactive  decay.  To  avoid  this  stiffness,  solutions  are 
most  commonly  found  by  making  approximations.  In  this  typical  implementation, 
isotopes  are  evaluated  differently  based  on  the  magnitude  of  the  half  life  in  relation  to 
other  isotopes  in  the  chain  and  the  position  of  the  isotope  in  the  chain  (the  first  isotope  in 
the  chain  versus  an  isotope  near  the  end  of  the  chain). 

For  some  applications,  models  find  gamma  activity  by  using  an  approximation 

_ i  o 

that  the  gamma  activity  decays  as  t  '  .  This  approximation  is  used  for  any  arbitrary 
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gamma  energy  range  as  well  as  total  activity.  The  new  model  seeks  to  improve  upon  this 
approximation  by  calculating  the  activity  directly  for  each  time  step  for  any  given  energy 
bin. 

This  chapter  outlines  the  motivation  for  developing  a  new  method  for  use  in 
radioactive  decay  calculations,  the  goal  of  the  research,  the  scope  of  the  problem  I 
attempted  to  solve,  and  the  assumptions  made  in  developing  this  new  implementation  of 
the  exponential  moments  function. 

I.  A :  Motivation 

A  more  robust  method  of  calculating  radioactive  decays  is  desired  to  increase  the 
precision  of  the  results.  A  method  which  treats  all  isotopes  the  same,  regardless  of 
position  in  the  chain  or  relative  magnitude,  will  treat  decay  chains  uniformly  and 
therefore  more  predictably  and  with  greater  precision. 

It  is  important  to  know  the  quantities  of  each  isotope  remaining  after  a  given  time 
more  precisely  than  available  with  current  methods.  This  will  allow  for  more  precise 
calculation  of  the  activity  of  the  radioactive  isotopes  and  lead  to  a  better  understanding  of 
high  energy  gammas  emitted  following  the  neutron  induced  fission  of  uranium  and 
plutonium  isotopes. 

The  Way-Wigner  approximation  is  that  the  dose  rate,  D ,  following  thermal 
fission  of  U-235  falls  off  as 

D(t)  =  D(l)r12  (1.1) 
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where  D(\)  is  the  dose  rate  at  any  unit  time[4:2].  We  examine  the  accuracy  of  this 
approximation  as  applied  to  various  fuels,  spectra  of  neutrons  that  induce  fission,  and 
gamma  radiation  energy  groups. 

I.  B :  Statement  of  the  Problem 

Precisely  model  the  time-dependent  gamma  spectrum  produced  by  the  neutron 
induced  fission  of  three  different  isotopes  with  three  different  energy  neutrons  while 
quantifying  the  gamma  activity  uncertainty  due  to  uncertainties  of  the  isotope  half  lives. 
Develop  computer  programs  that  are  useful  tools  for  performing  such  studies. 

I.  C:  Goal  of  the  Research 

The  goal  of  this  research  is  to  more  precisely  characterize  the  gamma  activity 
following  a  fission  event  by  using  the  exponential  moments  function  in  a  radioactive 
decay  code  to  calculate  the  quantity  of  isotopes  remaining  at  some  given  time  with 
precision.  This  research  also  attempts  to  quantify  the  propagation  of  uncertainty  in  the 
gamma  activity  caused  by  uncertainties  in  both  the  half  lives  of  isotopes  and  gamma 
decay  intensities. 

I.D:  Scope 

The  scope  of  this  research  is  limited  to  creating  a  decay  code  for  all  isotopes  with 
fewer  than  99  protons  (californium  and  below).  Isotopes  with  higher  Z  numbers  can  be 
used,  but  the  amount  of  information  available  for  elements  above  californium  drops  off 
significantly. 
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Fissionable  nuclei  considered  here  are  limited  to  uranium-235,  uranium-238,  and 
plutonium-239.  For  each  of  these  isotopes  three  different  incident  neutron  spectra  are 
evaluated:  thennal  neutrons,  fast  fission  spectrum  neutrons,  and  high-energy  neutrons. 
Because  U-238  does  not  fission  by  absorbing  a  thermal  neutron,  a  total  of  eight  different 
test  cases  are  evaluated. 

The  specified  time  step  for  the  code  can  be  any  number  between  1 .0  nanosecond 
and  25,000  years. 

I.E:  Assumptions 

The  fission  product  lists  were  taken  from  text  files  created  by  T.R.  England  and 
B.F.  Rider  and  made  available  on  the  Lawerence  Berkely  Laboratory  website  [5].  These 
files  contain  the  atomic  number  and  symbol  of  the  isotope,  as  well  as  an  indication  if  the 
isotope  is  in  its  ground  state  or  a  metastable  state.  For  each  metastable  state  listed,  I 
assumed  that  this  corresponded  to  the  most  excited  metastable  state  found  in  the  NuDat 
2.0  data  files.  If  an  isotope  listed  as  a  fission  product  does  not  appear  in  the  NuDat  2.0 
data,  I  assumed  that  it  instantly  underwent  successive  beta  decays  until  reaching  an 
isotope  for  which  there  is  data. 
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II.  Analysis  and  Approach 


A  functional  schematic  of  the  approach  used  to  solve  the  complex  decay  chain 
problem  is  shown  in  Figure  1 .  The  white  boxes  represent  the  formatted  text  files 
containing  user  input  data  and  code  output  files.  The  shaded  boxes  show  the 
macroscopic  processes  performed  by  the  code. 


Input  Code  Tasks  Output 


Figure  1  Functional  Code  Schematic 

The  rest  of  this  chapter  is  dedicated  to  a  discussion  of  radioactive  decay  theory 
and  examining  the  macroscopic  code  tasks  in  greater  detail. 

II.  A :  Special  Radioactive  Decay  Problem 

The  special  radioactive  decay  problem  describes  the  decay  of  an  initial  quantity  of 
a  single  radioactive  isotope  that  is  not  being  replenished  and  decays  through  a  single 
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chain  (without  branching).  Schematically,  this  can  be  represented  using  generic  isotopes 
1,  2,  3,  and  so  on  to  an  arbitrary  depth  ending  in  the  kth  isotope. 

1 — — — >2 — — — >3 — — — > - bht± — >k 

The  branching  fractions  (i.e.  b2  \  )  are  all  equal  to  unity  for  the  special  decay  problem. 

The  decay  of  the  initial  isotope  1  is  described  by  the  differential  equation 

(t)  +  AlNl  (t)  =  0  (2.1) 

where  \  and  Nx  (/)  are  the  decay  constant  and  quantity  of  isotope  1  at  time  t.  The 
differential  equation  for  the  quantity  of  isotope  2  is  similarly  given  by 

N1(t)  +  l1N1(t)  =  bu\Nx(t).  (2.2) 

This  method  can  be  generalized  for  the  decay  chain  of  arbitrary  length  k  as 

Nk(t)  +  \Nk(t)  =  bkjc-\A'k-\Nk-\(t)  (2.3) 


The  special  radioactive  decay  problem  for  the  sample  isotope  decay  can  be 
represented  as  a  lower-triangular  matrix  as  in  equation  (2.4): 


'  ^1(0  'I 

(  K 

0 

0  ••• 

0 

0  " 

'  Nx(t)  > 

N2(t) 

h. 

0  ••• 

0 

0 

N2(t) 

W) 

+ 

0 

^3 

0 

0 

N3(t) 

=  N(t)  +  AN(t)  =  0  .(2.4) 

4-t(0 

0 

1 

0 

0  0 

0 

Nk-i  (t) 

l  Nk(t)  J 

l  0 

0 

0  0 

h; 

{  Nk(t)  j 

One  solution  to  equation  (2.4)  is  of  the  fonn 

N(t)  =  e~AtN(0) .  (2.5) 

By  defining  the  transmutation  matrix,  or  T-matrix,  as 


T(t)  =  e~At, 


(2.6) 


6 


the  solution  to  equation  (2.4)  can  alternatively  be  expressed  as 

N(t)  =  T(t)N(  0).  (2.7) 

The  advantage  of  this  fonnulation  is  that  elements  of  the  T-matrix  are  calculated 
individually  and  the  solution  is  calculated  using  matrix  operations.  This  takes  advantage 
of  the  low  memory  requirements  for  individual  element  calculations  with  the  speed  of 
optimized  matrix  multiplications.  The  solution  method  used  to  calculate  the  values  of 
each  element  of  the  T-matrix  affects  the  efficiency,  accuracy,  and  precision  of  the  results. 

The  calculation  of  the  activity  is  straightforward  once  N(t)  is  known.  The 
activity  for  the  special  decay  problem  produced  by  any  isotope  is  defined  as  the  product 
of  the  decay  constant  and  the  quantity.  A  subset  of  the  total  activity  of  interest  in  this 
thesis  is  the  gamma  activity.  Using  the  concept  of  the  transmutation  matrix  introduced 
above,  the  gamma  activity  can  be  calculated  by 

A(t)  =  GT(t)N(0)  =  GN(t).  (2.8) 

Individual  elements  of  the  gamma  matrix  are  populated  by  calculating  the  product 
of  the  decay  constant  and  the  intensity  of  the  radiation.  The  radiation  intensity  indicates 
the  probability  of  a  specific  gamma  ray  being  emitted  for  a  given  decay  process.  Because 
the  radiation  intensity  is  given  for  each  gamma  radiation,  it  possible  to  sort  the  gamma 
activity  by  energy  bins  specified  in  the  code  by  the  user. 

II.  B:  Solution  Methods  of  the  Special  Problem 

I  investigated  four  different  approaches  to  solving  the  radioactive  decay 
differential  equations:  the  Bateman  solution  fonnula,  numerical  integration  of  the 
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ordinary  differential  equation,  matrix  exponentiation,  and  Mathew’s  solution  in 
exponential  moments  functions  [7]. 


1.  Bateman  Solution  Formula 

The  most  common  way  of  solving  the  differential  equations  that  describe 
radioactive  decay  is  known  as  the  Bateman  solution.  The  Bateman  solution  is  easily 
derived  using  Laplace  transforms  as  shown  in  numerous  nuclear  engineering  texts  (1:76- 
78).  The  Bateman  solution  formulates  the  solution  to  equation  (2.1)  as 

Nl(t)  =  Nl(  0)e~Alt.  (2.9) 

The  solution  to  the  second  differential  equation  begins  by  recognizing  that  the  right  hand 
side  of  equation  (2.2)  includes  /V,  (t)  which  is  known  explicitly  from  equation  (2.9). 

Making  this  substitution  and  algebraically  manipulating  the  result  leads  to  the  solution  for 
the  quantity  of  the  second  isotope 


(  -A, t 


N2(t)  -  b2  \\Ni{Q) 


~h.t 


(Aq  -/ij)  (2-J  ~Aq) 

Generically  applying  this  formula  to  the  last  (kth)  isotope  in  a  decay  chain  yields  the 
Bateman  solution  formula: 


(2.10) 


f  k-l 


Nk(t)  =  N1(  0) 


k 


~A:t 


nw.  s- 


V  1=1 


J 


i= 1 

i±j 


(1.11) 


The  principle  advantage  of  the  Bateman  solution  is  the  straightforwardness  of  the 
approach.  The  formulation  makes  it  easy  to  implement  into  a  computer  code  with  nested 
DO  loops  to  handle  the  summation  and  product  operations. 
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The  denominator  of  the  Bateman  solution  however  involves  a  subtraction  between 


the  decay  constants  of  the  isotopes  in  the  decay  chain.  When  the  decay  constants  are 
very  close  together,  catastrophic  cancellation  occurs  resulting  in  a  shrinking  denominator 
which  causes  the  loss  of  digits  of  precision.  Implemented  into  a  computer  code  directly, 
this  loss  of  precision  leads  to  significant  errors  in  the  final  result. 

2.  Numerical  Integration  of  Coupled  ODEs 

Another  approach  to  obtaining  the  answer  to  the  decay  problem  is  to  use  a 
numerical  integration  method  to  solve  the  coupled,  linear,  first  order  differential 
equations  directly.  For  small,  well  defined  problems,  the  primary  advantage  of  this 
solution  is  that  the  loss  of  precision  inherent  in  the  Bateman  equation  when  At  ~  A  •  is 

avoided  by  solving  the  differential  equations  directly. 

The  advantage  of  this  approach  is  offset  by  the  large  computational  resources 
which  are  necessary  to  simultaneously  solve  individual  decay  chains  and  the  stiffness 
which  results  when  At  and  AJ  differ  by  many  orders  of  magnitude.  Scaling  this  beyond 

the  scope  of  limited  test  problems  to  calculate  all  of  the  decays  of  interest  to  this  thesis 
would  be  prohibitive  given  the  3,443  unique  isotopes  and  more  than  68,000  unique  decay 
chains.  The  NDSolve  built-in  routine  in  Mathematica  [9]  was  used  to  numerically  solve 
the  differential  decay  equations  for  a  series  of  test  problems  representative  of  the  actual 
decay  isotope  data  using  up  to  200  digits  of  working  precision  and  with  decay  chains  up 
to  four  isotopes  deep.  Even  under  this  limited  test  condition  the  numerical  integration 
solution  displayed  numerical  instability  and  took  hundreds  of  times  longer  to  calculate 
than  other  methods. 
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3.  Matrix  Exponential  (Transmutation  Matrix)  Method 

The  matrix  exponential  method  is  used  to  solve  a  system  of  linear  first  order 
differential  equations  with  constant  coefficients.  This  method  represents  the  system  of 
equations  below 

N2(t)  =  b2^N,(t)-A2N2(t\  N2{  0)  =  0 

N;  it)  =  bk_u^N2  it)  -  4  Nk  (t),  Nk(  0)  =  0 
in  the  matrix  form 

N'(t)  =  CN(t),  N(0)  =  N°  (2.13) 

where  C  is  the  square  coefficient  matrix  of  order  k.  The  solution  to  equation  (2.13)  is  of 
the  form 

N{t)  =  etCN{  0)  (2.14) 


and  can  be  solved  through  use  of  a  matrix  of  the  change  basis  as  described  in  Matrix 
Theory  and  Linear  Algebra  [6:376-405].  ORIGEN  uses  this  type  of  formulation  to  make 
its  calculations  in  the  method  shown  in  equation  (2.15),  where  tC=  A  . 


.  a2  £> 

eAN  =  (I  +  A  +—+  —  +  ---)N 
2!  3! 

A  A  A 

=  N  +  AN  +  ^(AN)  +  j(j(AN))  +  --- 


(2.15) 


The  advantage  of  using  the  second  equation  is  that  the  computer  never  has  to  multiply 
anything  larger  than  a  matrix  on  a  vector,  reducing  the  computational  cost  of  this  method. 
The  series  can  be  expanded  to  as  many  terms  as  needed  to  achieve  the  desired  precision. 
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The  formulation  of  the  problem  as  a  transmutation  matrix  is  a  useful  construct.  In 
practice,  however,  perfonning  matrix  operations  on  a  square  matrix  of  order  equal  to  the 
total  number  of  isotopes  is  computationally  too  expensive  to  perform  on  a  desktop 
computer  without  making  the  approximations  used  in  ORIGEN. 


4.  Transmutation  Matrix  by  Exponential  Moments  Function 

The  solution  method  presented  in  this  section  was  provided  by  Mathews  [7].  It 
uses  that  exponential  moments  functions  that  were  developed  by  Mathews  and  his 
students  [8].  They  were  designed  in  the  context  of  developing  discrete  ordinates 
transport  methods  but  can  be  applied  to  the  problem  of  radioactive  decay  to  avoid  the 
catastrophic  loss  of  precision  of  the  Bateman  solution.  The  moments  function  in  its  most 
general  case  is  defined  as 

1  h  4-i 

Mn(Al,A2,...Ak)  =  \dh\dt2...  J  (2.16) 

0  0  0 

For  the  radioactive  decay  problem,  the  unforced  n  =  0  or  M0  solution  is  used. 

1  h  4-i 

MQ(\,l1,...?Lk)  =  \dtx\dt2...  J  dtke~A'he(^)t2...e(Ak-'~Ak)lk.  (2.17) 

0  0  0 

The  exponential  moments  function  solution  is  suited  to  this  application  because  it  can 
take  advantage  of  a  recurrence  relation  to  calculate  the  solution  for  a  chain  of  any  length. 

After  multiplying  both  sides  of  equation  (2.2)  by  e^1 ,  the  quantity  of  isotope  2  at 
any  given  time  can  be  determined  as  above  by  integrating  over  the  time  of  interest. 


[N^  +  ^N^t i)]^1  =bn\Nx{tx)e^  =| ~{N2{h)e^) 


(1.18) 
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(1.19) 


J  j-(N2{tx)e^  )dtx  =  N2(t2)e ^  -N2{ 0)  =  J  -^-b2^Nx(tx)e^dt- 


o 


o 


L 

dtx 


N2(t2)  ~  ^2,i\e  ^  2  J  1 


(1.20) 


But,  N\  (q)  is  known  from  equation  (2.9).  Substituting  this  back  into  equation  (1 .20) 


yields 


2  2 

N2(t2)  =  b2X\e~^h  J  dtxNx( 0)e~^he^1  =  b2xlxe~^2  J  dtxNx( 0)e“(i™?1  (1.21) 


The  equation  for  any  isotope  j  along  the  chain  follows  the  same  form  as  equation 


(1.21). 


W  =  hi,H;\He  ''jti  j  dtHNH(tH)e 


xjtj- 1 


(2.22) 


Substituting  solutions  for  Nj_x(tj_x)  into  equation  (2.22)  and  rearranging  terms  reveals  a 


nested  structure  to  the  solution 


0  A 


N j^j^  =  bj,j-ibj-x,j-2^j-2e~Xjt]  j dtj 


h-i 


J  dt j_2N j_2{t j_2)e 


VlO-2 


(2.23) 


that  can  be  extended  to  the  generic  case  of  k  isotopes 


fk- 1  h 


Lk- 1 


e~Zkt ytk_ie-^k-i-h)tk-i  J  dtk_2e~(Ak-2~Ak-l)tk-2...jdtxNA(0)e 
0  0  0 

(2.24) 


V»4  J 
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Multiplying  and  dividing  through  by  t  and  performing  a  change  of  variables  where 


ti-t  U  j 

J  sit  .  J 


r  ml  p 

,  =  tj  I,  and  J  —  =  duH  then  rearranging  terms  produees  equation  (2.25) 


k- 1  )  i  uk-\ 

Nk(t)  =  Nl(0)  Ylbi+lAt  e~Akt^duk_]e~Uk-'l~Akt)llk-'  J  duk_2e{^k- W4-2*-4'K-2  _ 


(2.25) 


Now  applying  the  exponential  moments  function  to  equation  (2.25)  in  the  solution  of  the 
differential  equations  describing  radioactive  decay,  the  quantity  of  the  kth  isotope  can  be 
represented  by 


Nk(t)  =  Nl(0)  nw*  e  **Mo[(4-i t-ht)i(\_2t-Xkt)i...\\t-Xkt)\.  (2.26) 


After  some  algebraic  manipulations  of  equations  (1.11)  and  (2.26),  the  Bateman  and 
exponential  moments  function  solutions  can  be  expressed  as  a  ratio  of  the  quantity  of 
isotope  k  to  the  initial  quantity  of  isotope  1 . 

N  (t\  (k- 1  k  -Af 

zA — 

7Vl(U^  4=1  )i= 


j=l  n(^-it) 

*■= i  7 

i*j 


(2.27) 


nvi^  e  4?mo [(4-fi - 40 > (4-2^-V)’-’ (4* -40] 


Equation  (2.27)  shows  that  the  exponential  moments  function  and  the  modified  Bateman 
equation  have  a  common  product  in  parentheses.  Dividing  through  by  the  common 
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product,  I  define  a  new  function,  F(/J) ,  dependent  upon  the  vector  of  decay  constants 
and  time: 


F(Xt)  = 


j=1  r i  -At) 

i= 1 


i*j 


(2.28) 


=  e  4?m0  [(4 - 4?) , (4_2t  - 40 - 40] 

Algebraic  manipulation  of  the  Bateman  solution  (Appendix  A)  shows  that  it,  and  hence 
F(Xt)  ,  is  invariant  with  respect  to  argument  order.  This  is  important  because  it  means 
that  the  decay  constants  can  be  sorted  by  size  before  performing  the  calculation.  By 
sorting  the  decay  constants  from  largest  to  smallest,  the  moments-function  form  of 
F(Xt)  can  be  rewritten  as 


F{Xt)  =  e  Smallest1  M 


^ excluding^  ^ smallest ^ 
smallest 


(2.29) 


This  ensures  that  every  argument  passed  to  the  exponential  moments  function  is  greater 


than  or  equal  to  zero.  If  the  arguments  passed  to  the  moments-function  are  all  non¬ 


negative,  then  0  <  M0  <  1 ,  eliminating  overflow  errors. 


The  final  fonn  of  the  equation  used  to  determine  the  quantity  of  isotopes  in  the 


decay  chain  at  time  t  is  therefore 


Ml 

NM 


f  k- 1  ^ 

g  ^smallest^ ^ 

nAi+i-v 

^excluding*  smallest^ 

V  i= 1  J 

smallest 

(2.30) 
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II.  C:  Generalization  to  the  Full  Problem 


To  this  point,  the  branching  fractions  have  been  assumed  to  be  unity  for  the 
special  radioactive  decay  problem.  In  reality,  many  radioactive  isotopes  have  multiple 
decay  mechanisms  resulting  in  branching  ratios  between  0  and  1 .  Figure  2  gives  an 
example  of  a  more  complex  decay  scheme.  Here  isotopes  1  and  2  have  multiple  decay 
mechanisms,  isotope  5  is  produced  from  more  than  one  parent  isotope,  and  isotopes  4  and 
6  are  stable. 


Figure  2  Complex  Radioactive  Decay  Chain 

The  system  of  differential  equations  describing  this  more  complex  decay  scheme  is 
shown  in  equation  (2.31). 


r  a 

0 

0  0 

0 

0 

\ 

rNx(t)' 

N2(f) 

~^2,\\ 

h. 

0  0 

0 

0 

N2(t) 

W) 

i 

0 

4  0 

0 

0 

N3(t) 

N4(t ) 

0 

-b42^2 

0  1 

0 

0 

N4(t ) 

W) 

0 

-^>5,2^2 

“^5,3^3  0 

^5 

0 

N5(t) 

Aw) 

l  0 

0 

0  0 

“^6,5^5 

1 

) 

wo. 

=  N(t)  +  AN(t)  =  0 


(2.31) 

It  follows  from  this  that  the  transmutation  matrix  must  then  be  equal  to 
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T(t)  =  e~At 


Ml  o  0 

^i(O) 

Ml  Ml  o 

w,(0)  N2(0) 

Ml  o  Ml 

Ni(0)  N3(  0) 

Ml  Ml  o 

A^O)  N2(0) 

Ml  Ml  Ml 

NM  N2(  0)  7V3(0) 

Ml  Ml  Ml 

N,(  0)  N2(0)  N3(  0) 


0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

Ml 

Ns(  0) 

Ml 

n5(  0) 


0 

0 

0 

0 

0 

1 


(2.32) 


where  the  individual  elements  are  calculated  using  equation  (2.30). 

The  generalization  of  the  problem  using  6  isotopes  in  the  example  above  is 
directly  scalable  to  the  full  scope  of  the  thesis  with  3,443  isotopes.  The  T-matrix  for  the 
full  problem  is  sparse  (composed  of  a  vast  majority  of  elements  with  a  value  of  0)  and  is 
not  a  lower  triangular  matrix. 

Once  the  T-matrix  is  created,  the  quantity  of  each  isotope  at  time  t  is  detennined 
using  equation  (2.7)  and  taking  advantage  of  the  MATMUL  built-in  function  in  the 

Fortran  compiler.  Following  the  detennination  of  N(t)  ,  the  calculation  of  the  gamma 
activity  at  time  t  is  straightforward  using  the  gamma  activity  matrix  introduced  for  the 
special  decay  problem  above.  However,  elements  of  G  are  now  calculated  as: 

r  \ 


(2.33) 


k  \  l  J 


gjj  is  the  gamma  activity  due  to  radioactive  decay  in  energy  bin  i  from  isotope  j.  Aj  is  the 


decay  constant  for  isotope  j  .  bk  ■  is  the  branching  ratio  for  the  decay  from  isotope  j  to 
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any  isotope  k  that  results  in  a  gamma  ray  emission.  Ij  k(Ei)is  the  intensity  of  the 
radiation  line  /  in  energy  bin  i  resulting  from  the  decay  of  isotope  j  into  isotope  k  . 
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III.  Implementation 


The  gamma  ray  activity  code  is  made  up  of  five  main  sub  tasks:  the  depth- first 
search  used  to  identify  all  of  the  daughter  isotopes  in  a  chain,  the  decay  algorithm 
implementing  the  exponential  moments  function,  the  generation  of  the  T-matrix,  the 
calculation  of  the  quantities  of  each  isotope,  and  finally  the  calculation  of  the  gamma  ray 
activity  rate.  The  rest  of  this  chapter  discusses  the  data  necessary  to  run  the  code  and 
then  examines  the  functioning  of  each  one  of  the  5  sub  tasks  in  detail. 

III.  A :  Input  Data 

The  data  required  to  run  the  code  are:  radioactive  decay  data  including  the  decay 
mode  and  half  life  with  uncertainties  for  every  isotope,  gamma  decay  data  including 
gamma  energies,  decay  modes,  and  intensities  with  uncertainties,  and  a  list  of  initial 

quantities  of  isotopes,  N(0)  .  The  decay  code  is  designed  to  work  with  any  set  of  input 
data  as  long  as  it  is  in  the  proper  fonnat.  The  individual  processes  were  verified  using 
test  problem  sets  containing  completely  generic  data,  as  will  be  discussed  in  more  detail 
in  Chapter  IV.  This  means  that  performance  of  the  code  is  solely  dependent  upon  the 
quality  of  the  data  that  is  used  to  characterize  the  problem. 

1.  Isotope  Decay  Information 

Isotopic  decay  data  was  taken  from  NuDat  2.0,  an  interactive  web  based  program 
developed  by  the  National  Nuclear  Data  Center  (NNDC).  NuDat  2.0  produces  a  text 
output  of  all  isotope  decays  with  data  obtained  from  the  Evaluated  Nuclear  Structure 
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Data  File  (ENSDF).  The  NNDC  maintains  the  ENSDF  file  and  currently  lists  3,165 
different  isotopes  and  more  than  80,000  gamma  rays.  The  difference  in  the  number  of 
isotopes  maintained  by  the  NNDC  and  the  number  of  isotopes  used  in  the  decay  code 
(3,443)  is  due  to  the  accounting  of  metastable  isotope  states.  The  NNDC  counts  all 
metastable  states  as  one  isotope,  whereas  the  decay  code  assigns  a  different  isotope 
number  to  each  metastable  state  for  tracking  purposes. 

An  example  of  the  NuDat  2.0  output  is  below  in  Table  1.  The  only  changes  made 
to  the  table  from  the  original  data  are  that  columns  containing  infonnation  on  the  spin 
and  uncertainties  in  the  mass  excess  energy  and  natural  abundance  have  been  omitted  for 
readability.  The  uncertainties  are  contained  in  the  text  version  of  the  half  life  and  given 
in  the  Nuclear  Data  Sheets  style  [5],  In  the  case  of  tritium,  the  half  life  is  12.32  ±  0.02 
years. 


Table  1  Sample  NuDat2  Output  for  Hydrogen,  Deuterium,  and  Tritium 


A 

Element 

Z 

N 

Mass  Exc 
(MeV) 

T  Vi  (txt) 

T  Vi  (s) 

Abundance 

(%) 

Decay 

Mode 

Branching 
Ratio  (%) 

1 

H 

1 

0 

7.289 

STABLE 

Infinity 

99.99 

2 

H 

1 

1 

13.136 

STABLE 

Infinity 

0.01 

3 

H 

1 

2 

14.95 

12.32  Y  2 

3.89E+08 

0 

B- 

100 

Several  problems  were  encountered  in  using  the  NuDat  2.0  output  data  for  a 
decay  code.  The  purpose  of  the  database  is  to  provide  as  much  infonnation  to  a  user  as 
possible,  while  for  the  purposes  of  writing  a  decay  code,  the  quantity  of  data  is  not  the 
most  important  factor.  What  instead  is  needed  is  the  infonnation  that  produces  complete 
decay  chains.  For  example,  Ni-48  undergoes  electron  capture  to  Co-48;  however  no 
isotope  information  is  available  on  Co-48.  The  best  way  I  found  for  dealing  with  this 
problem  was  to  create  a  dummy  isotope  with  zero  protons  and  zero  neutrons.  All  decays 
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for  which  no  daughter  isotope  information  exist  default  decay  into  this  imaginary  isotope. 
This  ensures  that  the  correct  amount  of  the  parent  isotope  decays  away  with  each  time 
step.  However,  this  treatment  of  unknown  daughter  isotopes  effectively  allows  leaking 
from  the  system  in  that  the  mass  of  recognized  isotopes  at  t  =  0  may  not  be  equal  to  the 
mass  of  the  recognized  isotopes  at  t  =  q  .  To  track  this  problem  an  error  log  is  created  for 

each  program  run  that  lists  the  isotopes  with  missing  decay  daughter  infonnation. 

Another  problem  with  the  data  is  the  half  life  uncertainties.  The  uncertainties  are 
given  as  a  a ,  or  standard  deviation,  for  a  Gaussian  distribution  [10].  However,  the 
uncertainties  listed  for  isotopes  can  not  be  a  cr  because  half  lives  can  not  be  less  than 
zero.  When  the  relative  uncertainty  is  small  compared  to  the  half  life,  this  approximation 
may  be  adequate  because  many  multiples  of  the  uncertainty  must  be  subtracted  from  the 
listed  half  live  to  create  a  negative  half  life  and  this  possibility  is  exceedingly  unlikely. 
However,  66  different  isotopes  have  half  lives  and  uncertainties  within  3crof  0.  Table  2 
lists  the  three  fission  product  isotopes  within  three  standard  deviations  of  0.  This 
problem  is  handled  using  rejection  sampling  and  is  discussed  in  greater  detail  in  Chapter 
V. 


Table  2  Fission  Product  Isotopes  with  Large  Relative  Uncertainties 


A 

Symbol 

Half  Life  (s) 

Uncertainty  (s) 

66 

Cr 

0.010 

0.006 

70 

Co 

0.50 

0.18 

76 

Ni 

0.24 

0.24 

2.  Gamma  Radiation  Data 

Gamma  radiation  infonnation  is  also  available  through  NuDat  2.0.  The  first  four 
data  fields  are  the  same  as  in  the  isotope  data  file.  Additional  fields  include  infonnation 
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on  the  decay,  radiation  type,  radiation  energy,  intensity,  and  dose,  each  listed  with  the 
appropriate  uncertainties.  A  sample  of  the  gamma  input  data  follows. 


Table  3  Sample  Input  Gamma  Data 


A 

Sym 

Z 

N 

Decay 

Mode 

T  !4  (txt) 

T  ‘A  (s) 

Rad 

Type 

RadE 

(keV) 

cr 

Rad 

Intensity 

(%) 

cr 

8 

He 

2 

6 

B- 

119.1  MS  12 

0.119 

G 

980 

84 

10 

7 

Be 

4 

3 

EC 

53.22  D  6 

5E+06 

G 

477.6 

20 

10.44 

4 

8 

B 

5 

3 

B+ 

770  MS  3 

0.77 

G 

511 

12 

12 

The  gamma  data  specifies  what  type  of  decay  produces  the  gamma  radiation 
(internal  transition,  fT ,  [i 1  ,  etc.).  The  code  treats  initial  quantities  of  metastable  states 
correctly  in  that  internal  transition  decays  and  the  associated  gamma  emission  are 
counted.  Subsequent  decay  mechanisms  in  the  code  do  not  populate  metastable  states 
because  all  decays  are  assumed  to  produce  a  nucleus  in  the  ground  state.  This  approach 
was  selected  because  of  the  lack  of  data  completely  describing  the  decay  process  through 
metastable  states.  While  this  introduces  errors  into  the  code  results,  most  metastable  state 
half  lives  are  very  short  compared  to  the  half  lives  of  the  ground  state  nucleii,  so  the  error 
is  not  intolerable. 


3.  A(0)Data 

Fission  product  yield  data  for  uranium-235,  uranium-238,  and  plutonium  239 
were  taken  from  work  by  T.R.  England  and  B.F.  Rider  [5],  They  define  the  thennal 
neutron  spectra  as  neutron  energies  which  are  in  thermal  equilibrium  with  the 
surroundings  present  in  a  light  water  reactor.  Fast  neutron  fission  products  are  the  yields 
caused  by  fission  spectrum  neutrons.  High  energy  neutrons  come  from  the 

\H  +  \h  — >  2 He  +  n  (D-T)  reaction  and  have  energy  of  14.07  MeV  [3:29]. 
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I  chose  to  use  fission  product  yields  as  population  sets  for  the  initial  values  of 


AfiF)  because  of  the  large  number  of  different  isotopes  produced  in  these  reactions.  The 
fission  of  U-235  by  thennal  neutrons  produces  776  different  isotopes.  This  provides 
complex  data  sets  with  well  documented  results. 

III.B:  Decay  Chain  Identification 

The  depth-first  search  is  initiated  by  sending  the  decay  chain  algorithm  a  starting 
isotope.  If  the  isotope  is  stable  then  the  algorithm  exits.  However,  in  the  case  of  a 
radioactive  isotope,  the  algorithm  identifies  the  number  of  mechanisms  by  which  the 
isotope  decays.  Each  one  of  the  decay  mechanisms  results  in  a  separate  recursive  call  to 
the  decay  chain  algorithm,  tracking  the  new  generation  information  as  well  as  adding  the 
decay  constant  and  the  branching  ratio  for  that  decay  to  global  storage  vectors.  The 
algorithm  recursively  calls  itself  until  a  stable  isotope  is  reached  at  the  end  of  the  chain. 

Upon  completing  the  chain  by  finding  a  stable  isotope,  the  vectors  containing  the 
branching  ratios  and  decay  constants  are  passed  to  the  decay  calculation  routine.  Once 
the  decay  calculations  are  perfonned,  the  algorithm  backs  up  one  level  and  either 
completes  the  process  for  the  next  branching  mechanism  or  passes  a  new  branching  ratio 
and  decay  constant  vectors  to  the  decay  calculation  algorithm  (now  each  reduced  in  size 
by  one). 

The  performance  of  the  depth  first  search  is  illustrated  by  returning  to  the  sample 
decay  scheme  introduced  in  Figure  2.  Figure  3  lists  each  decay  chain  for  the  sample 
complex  decay  scheme  in  the  order  in  which  they  are  identified.  The  figure  also  shows 
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the  branching  fraction  vector  and  decay  constant  vector  passed  to  the  decay  calculation 
subroutine. 
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Figure  3  Depth  First  Search  Chain  Identification 


III.  C:  Exponential  Moments  Function  Implementation 

Before  making  a  call  to  the  exponential  moments  function  routine,  the  At  vector 
which  consists  of  all  positive  and  decreasing  values,  must  be  created  by  multiplying  each 
decay  constant  in  the  chain  by  the  time  of  interest.  For  the  sake  of  clarity,  I  now  define  a 

new  vector  x  where 


x  ^ excluding ^  ^ 

smallest 


'smallest 


t . 


(3.1) 


The  x  vector  is  created  by  sorting  At  by  size,  and  subtracting  the  smallest  value  from 
every  other  value  in  the  vector.  Once  this  is  accomplished,  the  exponential  moments 
function  routine  is  called  and  the  solution  method  is  based  on  the  number  of  arguments  in 
x . 
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1.  One  Argument 


The  simplest  call  to  the  exponential  moments  function  solution  is  for  a  vector  of 
length  one  and  is  given  by 

Mq(x)  =  M0 (xj )  =  -  - — .  (3.2) 

*1 

Implementing  this  into  the  code  is  straightforward,  but  even  here  care  is  taken  to 
minimize  any  loss  of  precision.  Clearly,  if  x,  is  large  the  numerator  approaches  1.0 . 

Computational  time  can  be  saved  by  avoiding  the  expensive  exponential  calculation 
when  the  numerator  will  equal  1.0  to  the  specified  precision.  For  double  precision  (16 
digits  of  precision),  this  occurs  for  values  of  jq  >  36.8414 .  The  code  uses  the 
approximation  in  equation(3.3)  for  all  instances  where  x,  >  38.0  . 

M0(x1)  =  -  (3.3) 

*i 

At  the  other  extreme  is  values  of  x,  which  are  small.  Small  values  can  cause  a  loss  of 

precision  in  the  numerator  as  the  exponential  term  approaches  1 .0  which  is  amplified  by 
then  dividing  by  the  small  number.  To  avoid  this  a  series  expansion  is  perfonned  up  for 
values  of  Xj  <  0. 1 .  For  instances  where  x,  is  neither  small  nor  large,  the  default 
exponential  moments  function  formulation,  equation  (3.2),  is  used. 

2.  Two  Arguments 

For  a  vector  of  size  two,  the  moments  function  can  be  expressed  as  a  relationship 
between  the  exponential  moments  function  of  the  individual  arguments.  The  two 
arguments  are  sorted  such  that  xl  <=  x2 . 
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(3.4) 


Mq(x)  =  M()(x],x2)  = 


Mq(x2 )  MpjxO 


*1  ~x2 

The  same  logic  used  in  the  single  argument  case  above  can  be  extended  here,  and  if 
x2  >38.0  (since  x  is  sorted,  this  also  means  that  xx  >  38.0 )  then 


M0(x1,x2)  =  —  .  (3.5) 

*1*2 

Likewise,  if  both  arguments  are  small,  the  Taylor  series  expansion  is  used  to  cancel  terms 
and  avoid  loss  of  precision. 

There  is  also  a  new  case  that  must  be  investigated.  If  xl  and  x2  are  neither  both 

small  nor  both  large,  then  a  case  may  arise  where  they  are  close  together,  and  if  not  dealt 
with  correctly,  leaves  the  same  problem  which  is  so  destructive  to  the  precision  of  the 
Bateman  equation.  Here  though  the  exponential  moments  function  formulation  can  be 
re-expressed  using  a  change  of  variable  x2  =xx  +  dx  to  eliminate  the  subtraction  in  the 
denominator. 


M„(x„  x,  +dx)  =  M«(X 1  +A>-M°  <*■>  =  tfp(»l)-r"M0W 
011  xx-(xx+dx)  x2 

3.  Three  Arguments 

If  all  three  of  the  arguments  are  large,  the  inverse  product  approximation  is  used. 
If  all  three  of  the  arguments  are  small,  a  series  expansion  is  used.  For  the  case  of  three 
arguments  where  either  two  of  the  arguments  or  all  of  the  arguments  are  close  together, 
the  following  equation  is  used  to  evaluate  the  arguments  once  the  arguments  are  sorted 
into  increasing  order  ( xx  <=  x2  <=  x3  ). 
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M0 (JCJ , X, , x3 )  =  Mo(Xl’X2)~g 


(3.7) 


In  the  general  case,  the  3ld  kind  exponential  moments  function  is  handled  by  equation 
(3.8). 


i,  , . ,  (l-ifWo(h)-wo(h)  +  #o(h)  ..  x2  ~ xi 

ivi o fXj , x2 , x3;  -  -  7-  —  )/-' 


x(t-%)(x3-xi  y 


(3.8) 


X3-X1 


4.  N  Arguments 

For  all  times  of  interest  in  the  implementation  of  this  code,  there  are  never  more 
than  three  not  all  large,  not  all  small,  close  together  arguments.  Therefore,  given  a  chain 
of  n  decays  the  decay  calculation  can  be  expressed  as  a  recursive  use  of  the  exponential 
moments  function,  dividing  x  into  subvectors  of  the  original  arguments.  Assuming  that 
the  arguments  are  sort  such  that  xl  >=  x2  >=■■■  >=  xn , 


M0(x)  =  M0(x1,x2,...,x„) 


M0(x1,...,x„_1)-M0(x2,...,x„) 


(3.9) 


x„-x 1 

Similarly  to  the  specific  argument  cases  of  the  exponential  moments  function  introduced 
above,  when  the  arguments  are  either  all  large  or  all  small,  an  inverse  product  or  a  series 
expansion  is  used.  If  neither  of  these  conditions  are  met,  the  moments  function  is  then 
called  recursively  until  the  arguments  passed  into  the  moments  function  are  all  large,  all 
of  the  arguments  are  small,  or  vectors  of  length  one,  two,  or  three  are  the  only  size 
vectors  remaining. 
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III.D:  T -matrix  Generation 


The  purpose  behind  the  T-matrix  is  to  be  able  to  perfonn  a  matrix  multiplication 
on  a  vector  of  the  initial  quantities  of  isotopes  to  produce  the  remaining  quantities  of  all 
isotopes.  The  T-matrix  is  initialized  as  an  identity  matrix  with  order  k  (where  k  =  the 
total  number  of  isotopes).  For  a  stable  isotope  number  j  where  1-7-  k ,  the  T-matrix 
value  will  remain  1.0  at  position  (j,j)  in  the  matrix  and  no  calculations  are  performed  on 
any  element  in  the  isotope  column. 

Each  radioactive  parent  isotope  is  passed  to  the  depth  first  search  routine.  Once  a 
chain  is  completely  identified,  the  information  is  passed  to  a  decay  calculation  routine 
which  calls  the  exponential  moments  function.  The  value  along  the  diagonal  is  calculated 
directly  for  every  radioactive  isotope  using 


Njit)  _,t 

— - - =  e  1 

*y(0) 


(3.10) 


Values  in  the  T-matrix  in  general  will  always  be  between  0.0  and  1.0.  The  only 
exceptions  to  this  are  decay  processes  that  result  in  more  than  one  of  the  same  type  of 
particle.  For  example,  lithium-3  undergoes  a  decay  which  results  in  the  emission  of  two 
protons,  effectively  undergoing  a  spontaneous  fission  with  three  hydrogen  atoms 
produced  by  each  decay  of  a  lithium-3  atom.  In  this  case,  the  T-matrix  value  (hydrogen- 
1,  lithium-3)  would  approach  3.0  for  any  time  of  interest  much  greater  than  the  half  life 

of  lithium-3  (  7.56*1 0  23  s). 

Decay  mechanisms  that  result  in  the  emission  of  neutrons,  protons,  or  alpha 
particles  have  these  particles  tracked  as  well.  For  example,  decay  by  proton  emission  to  a 
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daughter  isotope  is  also  recorded  as  decaying  to  hydrogen.  This  ensures  that  the  baryon 
count  is  conserved  throughout  the  decay  code. 

III.E:  Calculating  N(t) 

Once  the  T-matrix  has  been  created  for  a  time  step  t,  finding  the  quantity  of  each 
isotope  at  time  t  is  straightforward.  First  introduced  in  equation  (2.7),  calculating  N(t )  is 
computed  directly  through  matrix  multiplication  of  the  T-matrix  and  the  initial  quantities 
vector,  ATO) . 

N(t)  =  T(t)N(  0)  (3.11) 

An  advantage  of  using  the  T-matrix  is  that  several  problems  can  be  investigated 

simultaneously  as  long  as  each  problem  is  specified  by  its  own  ATO)  vector.  The 

number  of  vectors  of  initial  quantities  of  atoms  is  represented  by  columns  of  data  in  an 

initial  quantity  matrix.  To  demonstrate  this  ability,  all  eight  problem  sets  I  wished  to 

investigate  were  contained  in  separate  columns  of  the  ATO)  matrix  as  in  Figure  4. 

N(:,  1)  =  ATO) :  U-235, thermal  neutrons 
N(:,2)  =  N(0) :  U-235, fast  neutrons 
N(:, 3)  =  ATO) :  U-235, high  energy  neutrons 
N(:,  4)  =  N(  0) :  U-238,fast  neutrons 
N(:,5)  =  ATO) :  U-238,high  energy  neutrons 
N(:,  6)  =  ATO) :  Pu-239,thennal  neutrons 
iV(:,  7)  =  ATO) :  Pu-239,fast  neutrons 

N(:,  8)  =  Ar(0) :  Pu-239,high  energy  neutrons 

Figure  4  Matrix  Storage  of  Eight  Fission  Product  Yield  Vectors  at  t=0 
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The  series  of  N(0)  vectors  can  be  combined  into  an  N(0)  matrix,  changing 
equation  (3.1 1)  to 

N(t)  =  T(t)N(0).  (3.12) 

The  Fortran  built-in  function  MATMUL,  when  called  with  the  arguments  of  the  T-matrix 
and  ATO)  matrix,  performs  this  calculation  more  efficiently  than  possible  through  direct 
coding  of  the  function. 

III.F:  Calculating  A(E,t ) 

Gamma  radiation  information  was  taken  from  the  NuDat  2  database  by  looking 
for  all  decays  that  produced  gammas.  The  80,125  results  produced  by  this  search  are 
used  to  populate  the  data  for  calculation  of  the  gamma  activity  produced  by  the  range  of 
isotopes  of  interest  to  this  thesis.  Gamma  energies  in  the  database  range  from  76.5  eV  to 
1 1.2589  MeV,  with  the  bulk  of  these  occurring  at  the  lower  energies.  For  example,  there 
are  only  three  gammas  of  over  10.0  MeV  in  the  database,  and  all  of  them  are  produced 

through  the  J3 1  decay  of  sodium-20.  The  distribution  of  gammas,  broken  down  into  200 
keV  bins  is  shown  in  Figure  5  and  Figure  6  below. 
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y  Energy  (keV) 

Figure  5  Histogram  of  Database  Gamma  Energies  (200  keV  bins) 
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y  Energy  (keV) 

Figure  6  Histogram  for  Gamma  Energies  >5  MeV  (200  keV  bins) 

Any  number  of  energy  bins  can  be  specified  and  defined  by  the  user  in  the  input 
file,  ‘ProblemData.txt’.  The  energy  bins  are  defined  by  a  minimum  and  maximum  bin 
energy  given  in  keV.  Applying  the  energy  bin  construct  to  the  matrix  fonnulation  of  the 
decay  problem  results  in  A (E,t)  =  G{E)T(t) N(0)  =  G(E)N(t)  . 


30 


IV  Verification  Process 


This  research  attempts  to  increase  the  precision  of  current  exponential  decay 
schemes  by  implementing  the  exponential  moments  function.  This  means  that  a  thorough 
verification  process  had  to  be  designed  and  implemented  to  ensure  that  the  new  code  was 
producing  the  high  precision  results  anticipated.  Verification  was  accomplished  using  a 
comprehensive  set  of  test  problems  that  covered  the  entire  range  of  potential  anticipated 
inputs  to  the  exponential  moments  function  as  well  as  various  decay  chains  and  gamma 
radiation  decay  energies  and  intensities.  Mathematica  was  used  as  a  benchmark  for  the 
solutions  by  performing  calculations  with  greater  than  machine  precision. 

IV.  A:  Defining  the  Range  of  Calculation  Inputs 

The  range  of  the  test  problems  used  for  verification  was  bounded  by  the  extremes 
of  the  decay  constants  found  in  the  nuclear  data  tables  and  by  the  anticipated  range  of 
time  steps.  Isotopes  of  interest  range  from  hydrogen  (Z=l)  to  californium  (Z=98).  The 

shortest  half-life  in  the  problem  scope  is  Li-3  with  a  half  life  of  7.56*  10“  s  .  At  the  other 
end  of  the  spectrum  are  the  stable  nuclides,  each  with  effective  half  lives  of  infinity.  In 
practical  terms,  the  exponential  moments  function  will  be  validated  using  the  longest 

lived  radioactive  isotope,  Ge-76,  with  a  half  life  of3.79*10  ~s  .  Therefore,  the  decay 

constant,  X ,  is  bounded  between  1. 83 *10"33s"1  and  9.1 7*  1021s-1 .  Tests  were  also 
performed  using  X  =  0  for  the  stable  isotopes,  but  values  between  the  0  and 

1.83*10“33s“1  were  not  explored. 
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Verification  testing  is  also  limited  to  time  steps  ranging  between  a  nanosecond 

and  25000  years  (7.94*10n  s).  Time  steps  greater  than  this  were  not  tested  but  are  not 
expected  to  produce  inaccurate  results.  This  is  because  the  longer  time  step  sends 
increasingly  more  arguments  to  the  large  argument  approximation  of  the  exponential 
moments  function  and  this  approximation  becomes  better  as  the  arguments  become 
larger. 

Because  the  exponential  moments  function  deals  with  the  differences  between  At 
values  as  in  equation  (3.1),  care  must  be  taken  to  ensure  that  not  only  the  range  of  the 
overall  isotopes  is  covered,  but  that  the  difference  between  the  values  is  also  addressed. 
Using  16  digits  of  precision,  any  difference  between  the  two  smallest  A  values  that 

exceeds  16  orders  of  magnitude  has  no  effect  on  the  values  of  the  argument  vector,  x  . 
This  sets  a  practical  limit  to  the  range  of  arguments  necessary  to  fully  test  the  calculation 
routines. 

The  input  values  can  now  be  specified  for  the  test  problems  needed  to  verify  the 
operation  of  the  exponential  moments  function  subroutines.  Values  that  populate  x  must 

have  values  that  range  between  1.83*10'42  and  2.89*  1044  and  cover  up  to  21  orders  of 
magnitude  (providing  a  buffer  beyond  the  16  digits  of  working  precision).  Additionally, 
because  the  exponential  moments  function  routine  handles  special  cases  for  argument 
vectors  of  size  of  one,  two,  and  three  and  a  general  case  to  break  down  larger  problems 
into  one  of  these  three  cases,  x ,  need  only  consist  of  up  to  three  values  to  completely 
exercise  all  of  the  calculation  paths. 

In  practice,  verifying  this  section  of  the  code  was  accomplished  by  specifying  a 
1.0  second  time  step  and  adjusting  the  decay  constants  to  the  necessary  values  to 
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compensate.  While  not  strictly  necessary  from  the  point  of  view  of  the  code,  it  greatly 
simplified  the  input  to  Mathematica  for  use  in  the  benchmark  calculations. 

IV.  B:  Creating  the  Calculation  Inputs 

The  input  for  the  decay  constants  used  in  the  decay  calculation  routine  had  to  be 
rational  numbers  so  that  Mathematica  could  perfonn  the  calculations  with  an  arbitrary 
precision  to  produce  the  benchmark  values.  To  satisfy  the  rational  number  requirements 
of  both  applications,  I  decided  to  define  all  decay  constant  by  an  integer  value.  Each 
integer  value  was  then  used  as  part  of  an  exponent  of  10  to  generate  the  necessary  values 
specified  previously.  The  four  decay  constant  values  were  calculated  using  the  following 
equations. 

h 

=  101001 

(4.1) 

h 

4  =101001 
/i4  =  1 0‘4 

IV.  C:  Verifying  the  Calculations 

Mathematica  was  chosen  to  benchmark  the  results  of  the  exponential  moments 
function  routine  because  it  allows  a  user  to  specify  the  number  of  working  digits  of 
precision  used  in  calculating  the  results.  This  increased  working  precision  comes  at  the 
price  of  increasing  the  calculation  time  significantly.  The  Mathematica  solution  was 
generated  by  passing  in  the  decay  constant  values  and  calculating  the  answer  using  the 
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Bateman  solution.  The  Mathematica  notebook  used  for  verification  is  attached  in 
Appendix  C. 

A  symmetric  relative  difference  (SRD)  was  used  to  compare  the  results  from  the 
two  methods.  The  SRD  is  defined  as 


$  SRD  Oh  z)  -  { 


0,  y  =  z  =  0 

It-zI 


(\y\+\z\)/2' 


else 


(4.2) 


where  y  is  the  result  produced  by  the  exponential  moments  function  routine  and  z  is  the 
result  produced  by  Mathematica.  I  chose  to  use  the  SRD  as  a  measure  of  perfonnance 
because  it  is  easy  to  interpret,  with  a  value  of  2.0  being  the  worst  case  scenario  and 
indicating  that  the  results  are  uncorrelated  between  the  two  methods.  Conversely,  an 
SRD  close  to  zero  indicates  close  agreement  between  the  two  methods. 

I  experimented  with  the  working  precision  of  Mathematica  during  the 
computation,  looking  for  the  fastest  computation  time  that  still  produced  the  smallest 
SRD.  I  found  that  only  after  increasing  the  working  precision  to  at  least  100  digits  did 

the  two  computational  methods  converge,  producing  SSRD  =  1 .0 147*1 0  1 1 .  The 

SSRd  sets  a  limit  on  how  much  disagreement  there  is  between  the  benchmark  and  the 

code  calculations.  The  value  indicates  that  the  answer  produced  by  the  code  is  good 
through  at  least  ten  digits.  Because  no  half  lives  are  known  with  uncertainties  to  more 
than  five  digits,  the  Monte  Carlo  propagation  of  uncertainty  reflects  only  the  uncertainties 
in  the  half  lives  and  not  uncertainty  due  to  calculations  in  the  code. 

This  worst  case  scenario  was  recorded  for  the  following  argument  set 
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X  = 


^0.147996^ 

0.147837 

0.147679 


(4.3) 


Three  close  together  not  too  small  {xsmaUest  >  0.1 )  arguments  produced  the  worst  case 


conditions  for  the  code. 


The  requirement  that  Mathematica  have  at  least  100  working  digits  of  precision 
implies  that  current  methods  for  computing  the  decay  of  radioactive  isotopes  using  the 
Bateman  solution,  even  if  performed  using  quadruple  precision,  may  not  be  accurate  to 
the  same  precision  as  the  double  precision  exponential  moments  function  solution. 


IV.  D:  Verifying  the  Depth  First  Search  Routine 

The  decay  chain  algorithm  uses  a  recursive  subroutine  to  perform  a  depth-first 
search  of  all  daughter  products  for  a  given  isotope.  Testing  this  algorithm  meant  creating 
test  problems  made  up  of  test  isotopes  with  the  necessary  decay  chain  complexity  to 
ensure  that  any  chain  that  was  present  in  the  actual  isotope  list  could  be  completely 
described. 

The  first  and  most  simple  problem  is  the  decay  from  isotope  “A”  to  isotope  “B”, 

A  —>  B  ,  where  A  is  radioactive  and  can  only  decay  to  isotope  B  which  is  stable.  From 
here,  the  sample  chains  grow  longer  (up  to  18  isotopes  in  length)  and  more  complex  with 
an  increasing  number  of  branches  (the  maximum  number  of  branches  tested  is  three).  A 
complete  list  of  the  test  decay  chains  is  in  Appendix  B. 

Once  each  isotope  list  is  read  in,  the  depth-first  search  is  called  and  the  decay 
chains  are  found.  After  a  chain  is  completely  identified  the  decay  calculation  subroutine 
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is  called.  The  decay  constants  passed  to  the  calculation  subroutine  are  then  compared 
against  a  list  of  the  expected  decay  constants  created  by  hand  for  each  test  problem.  By 
this  method,  every  time  a  decay  chain  calculation  is  called,  the  information  is  compared 
to  ensure  that  it  is  complete  and  accurate.  Any  errors  in  the  depth  first  search  are  output 
to  the  screen  so  that  the  user  can  evaluate  the  problem. 

The  depth  first  search  also  provides  an  opportunity  for  a  further  test  of  the 
exponential  moments  function,  and  the  results  calculated  during  the  depth  first  search 
verification  are  compared  again  against  the  Mathematica  solution  using  the  SRD. 
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V  Monte  Carlo  Estimation  of  Uncertainty 


A  Monte  Carlo  simulation  was  used  to  quantify  the  effect  of  half  life  uncertainties 
on  the  gamma  activity  produced  by  neutron  induced  fission.  This  was  done  by  using  an 
approximation  of  a  Gaussian  distribution  of  pseudo-random  numbers  to  perturb  the  half 
life  of  each  isotope  using  the  uncertainty  and  then  tracking  the  changes  in  the  quantity  of 
each  isotope  and  the  gamma  radiation  activity  as  a  result.  By  doing  this  1000  times,  the 
effects  of  the  half  life  uncertainties  can  begin  to  be  quantified. 

Because  of  the  large  relative  uncertainties  in  the  half  lives  of  the  isotopes,  the 
distribution  cannot  in  fact  be  Gaussian.  A  truly  Gaussian  distribution  requires  an  infinite 
number  of  standard  deviations  from  the  mean  to  cover  all  possible  outcomes.  For  real 
world  data,  a  Gaussian  approximation  may  be  appropriate  if  a  reasonable  ( 3cr ) 
distribution  of  the  data  encompasses  nearly  every  possibility.  However,  with  the  data 
provided,  some  isotopes  have  standard  deviations  which  are  large  compared  to  the  half 
life  (see  Table  1)  and  predict  negative  half  lives  within  the  3cr  tolerance. 

Recognizing  that  calling  the  half  life  uncertainty  a  standard  deviation  does  not 
adequately  describe  the  data,  I  chose  to  treat  the  distribution  as  Gaussian  and  use 
rejection  and  re-sampling  to  ensure  that  every  half  life  has  a  positive  value.  Clearly  this 
approach  only  provides  a  first  cut  at  quantifying  the  uncertainty  in  the  activity 
calculations,  but  this  should  still  be  sufficient  to  gain  a  better  understanding  of  how  the 
uncertainties  are  propagated  through  decay  chains. 

An  approximation  to  a  nonnal  distribution  of  random  numbers  was  created  by 
taking  the  sum  of  twelve  pseudo-random  numbers  uniformly  distributed  between  0  and  1 
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and  subtracting  6.0.  This  method  provides  a  good  approximation  to  a  nonnal  distribution 
near  the  center,  but  loses  fidelity  at  the  tails  of  the  distribution.  For  the  purposes  of  this 
thesis,  the  behavior  at  the  middle  of  the  distribution  was  the  important  component  and 
this  method  models  the  nonnal  distribution  well  enough  for  tracking  the  propagation  of 
uncertainty  in  the  isotope  data. 

Once  the  random  number  is  produced,  the  value  was  then  multiplied  by  the 
uncertainty  in  the  half  life  and  added  to  or  subtracted  from  the  half  life.  If  a  non-positive 
value  for  the  half  life  is  obtained  by  this  method,  another  random  variable  is  chosen  until 
a  positive  half  life  is  found.  Random  variables  were  chosen  using  the  built  in  random 
number  generator  in  Compaq  Visual  Fortran[2]. 
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VI.  Performance ,  Results ,  and  Analysis 


The  implemented  exponential  moments  function  code  combined  with  the  use  of  a 
transmutation  matrix  is  faster  than  the  differential  equation  solver  or  the  Bateman 
solution  with  the  specified  digits  of  working  precision  in  Mathematica.  Moreover,  the 
code  produces  answers  that  are  more  precise  than  other  methods  using  the  same  digits  of 
precision.  The  rest  of  this  chapter  will  outline  the  perfonnance  of  the  code  as  well  as 
show  results  from  the  fission  product  test  problems. 

VI. A:  Performance 

The  time  it  takes  to  produce  N(t)  and  A(t)  with  the  code  is  highly  dependent 

upon  the  time  of  interest.  The  overhead  of  reading  in  the  decay  data  to  populate  N(0) 
values  and  gamma  data  must  be  performed  only  once  with  each  data  set.  Then  the 
populated  data  arrays  are  stored  to  the  hard  drive  in  binary  format  for  quick  processing 
during  future  runs. 

The  major  factor  in  the  run  time  of  the  code  is  the  production  of  the  T-matrices. 
Once  a  T-matrix  is  generated,  it  is  saved  to  the  hard  disk  in  a  sparse  format  for  future  use 
which  greatly  speeds  any  subsequent  calculations  using  the  same  time  of  interest.  The 
time  required  to  create  a  T-matrix  from  scratch  depends  on  the  time  of  interest. 

At  very  short  times,  nearly  all  arguments  passed  to  the  Fortran  routine  are  small. 
When  this  is  the  case,  series  expansions  are  used  which  are  calculated  quickly  in  Fortran. 
For  very  long  times,  arguments  to  the  exponential  moments  function  all  become  larger 
than  38.0  and  the  T-matrix  is  generated  quickly  at  these  times  as  well.  However,  for  time 
steps  where  many  of  the  arguments  passed  to  the  exponential  moments  function  are 
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neither  all  large  nor  all  small,  calculating  the  results  becomes  computationally  costly  as 
the  number  of  EXP()  calls  reaches  a  maximum. 


T-matrix  Time  Step  (10ns) 

Figure  7  Time  to  Generate  T-matrix  for  Given  Time  Steps 
The  computer  used  to  execute  the  code  had  a  2.8  GHz  Pentium  D  processor  (a 
dual  core  chip  running  the  code  on  only  one  of  the  processors).  However,  for  longer  time 
steps,  it  becomes  more  advantageous  for  the  user  to  use  stored  T-matrices  to  reduce  the 
run  time.  For  this  reason,  a  newly  generated  T-matrix  is  automatically  saved  in  a  sparse 
format  once  it  is  created. 

These  results  are  indicative  of  the  computational  cost  for  the  current  algorithm 
which  is  only  optimized  to  preserve  precision.  Optimizing  the  exponential  moments 
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function  module  for  speed  with  long  argument  vectors  should  result  in  significant 
reductions  in  computer  run  time. 


VI. B:  A(E,t) 

I  arbitrarily  used  1 1  different  gamma  energy  bins  to  sort  the  output  of  the  decay 
code.  One  energy  bin  covered  all  gamma  energies,  while  the  remaining  ten  each  covered 
a  1  MeV  span  without  overlap  between  1  MeV  and  1 1  MeV.  The  gamma  activity  varied 
greatly  from  the  freshly  irradiated  fuel  state  1.0  nanosecond  after  irradiation  to  the  last 
calculated  time  step  of  just  over  25,000  years. 

The  plots  of  total  activity  show  the  overall  effect  of  the  combination  of  a  fuel  and 
a  specific  neutron  energy  spectrum.  As  is  shown  in  Figure  8,  the  total  activity  is 
dominated  by  the  lower  energy  gammas.  Therefore,  the  activity  at  higher  energy  bins 
must  be  examined  separately  to  completely  understand  the  decay  process. 
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Figure  8  Energy  Bin  Contributions  to  Activity  for  U-235  Thermal 


42 


21.0  ki.0 k)t0 1-60 1'  80  U  z0 1-  ©0 1-  s0  U  t^O  U  c0  U  2O  U  u0  U  oO  U  i.-0 1- 2-0  l-e-0 l-r-0  l-s-0 1- 9-O U z-0  l-s-O  U 6-0 1 


Because  of  the  large  ranges  of  time  and  activity  involved,  the  most  useful  way  of 
looking  at  this  data  is  to  break  it  into  two  separate  sections:  an  early  time  which  starts  at  a 
nanosecond  and  ends  around  10  seconds,  and  a  late  time  for  all  times  after  10  seconds. 
Also,  by  changing  the  axis  scaling,  different  features  of  the  activity  curves  are 
highlighted  and  examined  in  greater  detail. 

1.  Total  Activity 

At  times  shortly  after  irradiation  the  type  of  fuel  profoundly  affects  the  overall 
activity.  U-238  shows  significantly  more  activity  than  U-235,  which  is  more  active  than 
Pu-239.  The  activity  from  the  fission  of  Pu-239  with  a  thermal  neutron  spectrum  so 
closely  matches  the  activity  of  Pu-239  exposed  to  the  fast  fission  spectrum  that  the  lines 
are  nearly  indistinguishable  in  Figure  9. 
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Figure  9  Total  Activity,  1  nanosecond  to  10  seconds,  Log-Linear 
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As  shown  above,  the  activity  holds  nearly  constant  until  around  one  tenth  of  a 
second.  The  flatness  of  the  curves  is  due  to  only  very  minor  changes  in  activity  at  these 
time  steps.  This  indicates  that  the  process  is  dominated  by  either  isotopes  that  have  half 
lives  on  the  order  of  a  tenth  of  a  second  up  to  a  few  seconds  or  the  isotopes  are  in  secular 
equilibrium  with  a  parent  isotope  that  has  a  half  life  in  that  same  range  .  Figure  10  shows 
this  finding  more  directly  by  using  a  semi  log  plot.  On  the  semi-log  axis  configuration, 
the  plot  of  activity  of  a  single  isotope  is  represented  as  a  straight  line  proportional  in 
slope  to  the  decay  constant  of  the  isotope.  The  curves  shown  in  Figure  10  are  not 
perfectly  straight  because  the  activities  are  made  up  of  many  isotopes.  Each  of  the  top 
five  contributors  to  the  total  activity  is  from  isotopes  with  half  lives  between  0.2  seconds 
(Rb-96)  and  2.1  seconds  (Zr-99),  as  is  shown  in 
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Table  4. 
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Figure  10  Total  Activity,  1  nanosecond  to  10  seconds,  Semi-Log 
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Table  4  Major  Contributors  to  Total  Activity  at  4.64*10  5s  and  Half  Lives 


Fuel 

Incident 

Neutron 

Energy 

Primary 

Contributors 

Half  Life 

U-235 

Thermal 

Sr-96 

1.07  s 

Y-98 

0.548  s 

Sr-97 

0.429  s 

Fast 

Sr-96 

1.07  s 

Y-98 

0.548  s 

High 

Y-98 

0.548  s 

Xe-134m 

0.290  s 

Sr-96 

1.07  s 

U-238 

Fast 

Rb-96 

0.2028  s 

Y-100 

0.735  s 

Sr-96 

1.07  s 

Sr-98 

0.653  s 

High 

Sr-96 

1.07  s 

Y-98 

0.548  s 

Rb-96 

0.2028  s 

Y-100 

0.735  s 

Pu-239 

Fast 

Sr-96 

1.07  s 

Y-98 

0.548  s 

Zr-99 

2.1  s 

High 

Xe-134m 

0.290  s 

Y-98 

0.548  s 
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Figure  1 1  shows  a  semi-log  plot  of  the  total  time  to  show  that  the  activity  curves 
at  late  times  are  each  very  similar  in  slope.  The  flatness  of  the  activity  curves  at  the  end 
of  the  calculation  range  is  due  to  the  same  isotope  in  each  fuel,  Sb-126.  The  half  life  of 
Sb-126  is  only  12.35  days.  Here  Sb-126  is  in  secular  equilibrium  with  Sn-126  which  has 
a  half  life  of  230,000  years.  The  activity  due  to  the  decay  of  Sb-126  accounts  for  more 
than  90%  of  the  total  gamma  activity  in  each  decay  chain.  The  differences  in  the  activity 
between  isotopes  reflect  the  initial  allocation  of  fission  fragments  by  mass  number.  High 
energy  neutrons  create  fission  products  with  a  mass  number  of  126  with  more  frequency 
than  fissions  induced  by  neutrons  of  lower  energy.  The  correlation  of  this  relationship 
between  the  production  of  fission  products  with  atomic  numbers  of  126  and  the  activity  at 
25,000  years  is  0.992. 


49 


in 

CD 

Is- 

00 

CJ> 

T— 

T— 

T— 

T— 

T— 

T— 

T— 

<u 

CD 

CD 

CD 

CD 

CD 

CD 

CD 

Cl) 

CD 

CD 

CD 

T — 

T — 

T — 

T — 

T — 

T — 

T — 

T — 

T — 

T — 

T — 

T — 

(uojssjj/s/sAeoep  A)  AjjAjpv 


Figure  11  Total  Activity,  1  nanosecond  to  25,000  years,  Semi-Log 
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(uojssjj/s/sAeoep  A)  AjjAjpv 

Figure  12  Total  Activity,  1  ns  to  25000  years,  Log-Log 
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Figure  12  uses  logarithmic  axes  to  show  detail  in  the  activity  curves  over  the  time 
of  interest.  Here  we  see  that  the  activities  are  due  to  a  combination  of  fuel  and  incident 
neutron  energy.  For  each  fuel,  the  activity  is  highest  for  the  D-T  neutrons  by  several 
orders  of  magnitude.  For  the  same  incident  energy  neutrons,  the  fuels  are  Pu-239,  U-235, 
and  U-238  in  order  of  decreasing  activity.  The  areas  of  the  curve  that  are  largely  flat  are 
not  completely  straight  due  to  many  small  contributions  by  isotopes  with  half  lives  of  the 
same  order  of  magnitude  as  the  time  in  question.  The  portions  of  the  curves  with  large 
slopes  indicate  times  where  the  magnitude  of  the  half  lives  of  the  principal  contributors  to 
the  activity  are  of  the  same  order  of  magnitude  as  the  time.  For  example,  at  the  one  year 
point,  the  principal  contributor  to  the  total  activity  for  plutonium  and  the  high  energy 
neutron  fission  of  U-235  is  Tb-158  with  a  half  life  of  180  years.  The  remaining  fuels  (U- 
235  thermal,  U-235  fast,  U-238  fast,  and  U-238  high)  have  nearly  the  same  change  in 
activity  at  the  one  year  point  due  to  the  same  two  primary  contributors,  Nb-95  (half  life 
of  34.991  days)  and  Zr-95  (64.0  day  half  life),  which  together  account  for  more  than 
75%  of  the  total  activity  in  these  fuels.  This  explains  why  four  of  the  curves  are  basically 
flat  and  four  have  large  slopes  at  1.0  year  in  the  graph  above. 

Figure  13  shows  plots  of  the  gamma  power  (the  product  of  activity  and  the 
gamma  energy)  with  respect  to  time  for  each  fuel  and  neutron  spectra.  Additionally,  the 

_ i  o 

plot  overlays  the  Way-Wigner  approximation  oft'.  This  clearly  highlights  the  poor 
approximation  of  the  actual  dose  rate,  especially  for  the  D-T  neutron  fission  of  Pu-239. 
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Figure  13  Gamma  Power  and  the  Way-Wigner  Approximation 
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2.  1-2  MeV 


The  plot  of  short-time  gamma  activity  in  Figure  14  closely  resembles  the  shape  of 
the  total  activity  curves  in  Figure  9.  The  1-2  MeV  gammas  account  for  approximately 
only  one  fourth  of  the  total  activity.  As  in  Figure  9,  Figure  14  shows  that  the  early  time 
activities  are  more  dependent  on  fuel  than  neutron  energy,  although  the  Pu-239  high 
energy  fission  curve  is  now  nearly  indistinguishable  from  the  U-235  thennal  curve. 

Once  again,  the  list  of  major  contributors  to  the  1-2  MeV  activity  bin  shows  that 
all  have  half  lives  in  the  several  tenths  of  a  second  to  a  few  seconds  which  corresponds 
well  to  the  roll  over  of  all  of  the  curves  that  occurs  somewhere  after  0. 1  seconds.  The 
nearly  straight  curves  in  the  semi-log  plot  shown  in  Figure  15  shows  that  the  activities  of 
each  of  the  fuels  is  caused  by  multiple  isotopes,  each  with  similar  half  lives.  The  near¬ 
parallel  nature  of  the  curves  indicates  that  the  mixes  of  isotopes  from  each  fuel  sample 
are  likely  to  have  many  of  the  same  daughters. 

Figure  16  shows  the  entire  range  of  calculated  data.  Once  again,  this  plot  closes 
matches  the  plot  of  the  total  gamma  activity.  The  steep  slope  at  early  times  shows  nearly 

parallel  lines.  The  activity  in  this  bin  at  2.5  *  1010  seconds  (approximately  800  years)  is 
primarily  due  to  the  decay  of  Tb-158  (180  year  half  life).  Activity  at  the  end  of  the  time 
of  interest  is  dominated  by  the  decay  of  Sb-126,  just  as  was  the  case  for  the  total  activity. 
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Figure  14  Activity,  1-2  MeV  Gammas,  1  ns  to  10  s,  Log-Linear 


55 


Time  (s) 


_ 

CD 

03 

E 

E 

4—* 

.c 

, _ , 

.cz 

-CZ 

CD 

C/3 

03 

a; 

C/3 

o> 

C/3 

O) 

.C 

CD 

.c 

0! 

03 

4— 

03 

4— 

03 

_£Z 

03 

ID 

ID 

ID 

00 

oo 

CO 

CO 

CO 

CO 

CM 

i 

CD 

CM 

i 

CO 

CM 

i 

CO 

CM 

i 

CO 

CM 

CM 

CM 

i 

CM 

O 

3 

O 

D 

D 

D 

D 

D 

CL 

CL 

CL 

00 


CD 


CM 


O 


O 

O 


(uoissjj/s/sAeoep  A)  A}|A|pv 

Figure  15  Activity,  1-2  MeV  Gammas,  1  ns  to  10  s,  Semi-Log 
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Figure  16  Activity,  1-2  MeV,  1  ns  to  25,000  y,  Semi-Log 
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3.  2-3  MeV 


During  times  up  to  10  seconds  after  irradiation  of  the  fuel  source,  the  activity 
curves  once  again  closely  follow  the  shape  of  the  1-2  MeV  activity,  as  well  as  the  total 
activity.  However,  as  the  bins  continue  to  increase  in  energy,  there  is  a  decrease  in  the 
magnitude  of  the  activity  in  each  bin.  Whereas  the  1-2  MeV  accounted  for 
approximately  one  fourth  of  the  total  activity,  the  2-3  MeV  gammas  only  account  for 
around  a  tenth  of  the  total  activity.  The  curves  shown  in  Figure  17  show  the  2-3  MeV 
gamma  activity  curves.  The  graph  shows  that  the  difference  in  fuel  source  is  even  more 
pronounced  in  this  energy  bin.  Sr-97  and  Y-98  are  the  significant  contributors  for  each 
fuel  and  fission  neutron  energy  in  the  2-3  MeV  energy  bin. 
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Figure  17  Gamma  Activity,  2-3  MeV,  1  ns  to  10  s,  Log-Linear 
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Figure  18  Gamma  Activity,  2-3  MeV,  10  s  to  25000  y,  Log-Log 
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Figure  18  shows  the  gamma  activity  in  the  2-3  MeV  energy  bin  for  all  times  of 
interest.  This  is  the  first  gamma  energy  bin  for  which  the  activity  drops  below  an 

intensity  of  10  gamma  decays  per  second  per  fission  before  the  simulation  was 
completed.  I  have  set  this  as  a  lower  threshold  for  all  graphs  to  facilitate  comparisons 
among  the  plots.  In  all  cases,  the  slope  is  nearly  vertical  when  the  curve  crosses  this 
threshold,  and  there  are  no  defining  characteristics  below  this  threshold. 

Figure  18  also  shows  that  activities  for  all  fuel  sources  and  incident  neutron 
energies  are  nearly  indistinguishable,  with  the  exception  of  the  sub  1 .0  second  activities 

shown  in  Figure  17  and  a  separation  in  the  activities  between  105  and  106  seconds.  This 
second  area  is  magnified  in  Figure  19. 
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Figure  19  Gamma  Activity  with  Uncertainty,  2-3  MeV,  1  day  to  3  weeks 
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The  uncertainty  bars  in  Figure  19  are  standard  deviations  for  a  normal  distribution 
of  the  activity  over  1000  test  runs.  The  activity  results  were  each  checked  for  nonnality 
using  the  Kohnogorov-Smimov  test.  The  standard  deviations  show  that  the  fuels  are 
largely  distinguishable  within  lcr  after  approximately  two  days  of  decay  time,  but  that 
after  approximately  three  weeks,  all  of  the  curves  are  statistically  indistinguishable. 

The  activity  curves  for  U-235  with  fast  and  thennal  neutrons  are  nearly 
indistinguishable,  however  the  rest  of  the  fuels  and  neutron  sources  show  up  to  almost 

two  orders  of  magnitude  difference.  At  104  seconds  (2.8  hours),  different  sources  each 
have  the  same  major  contributors  to  the  total  activity,  albeit  in  varying  intensities.  The 
main  isotopes  at  this  time  are:  Kr-88  ( thay  =  2.84  days),  Cs-138  ( thay  =  33.41  minutes), 

La- 142(  thay  =91.1  minutes).  Each  of  these  isotopes  is  effectively  gone  by  one  week  and 

the  activity  becomes  dominated  by  isotopes  from  new  mass  chains.  At  2  *  105  seconds 
(2.3  days),  the  primary  contributors  to  the  activity  become  dominated  by  the  metastable 
state  of  Te-131  (thay  =  1.35  days),  1-132  (thay  =  2.30hours),  and  La-140 

(thaif  =1.68  days).  The  1-132  activity  is  produced  following  the  decay  of  Te-132 

( thay  =  3.26  days).  After  3  weeks,  the  activities  from  each  graph  join  again.  At  this 

point,  the  transition  is  complete  and  more  than  90%  of  the  activity  is  caused  by  La- 140  in 
each  isotope.  During  the  period  of  interest  above,  the  La- 140  decay  was  largely  the  result 
of  the  isotope  being  produced  directly  from  fission.  At  later  times,  it  is  in  secular 
equilibrium  with  its  parent,  Ba-140  ( thay  =  12.75  days). 
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4.  3-4  MeV 


Activities  in  the  3-4  MeV  energy  bin  only  contribute  about  2%  of  the  total  activity 
at  early  times  as  shown  in  Figure  20.  The  separate  fuels  are  once  again  readily 
identifiable.  At  1.0  second  into  the  decay,  the  activity  is  dominated  by  Y-97  for  all  test 
cases. 

—9  5 

Figure  21  depicts  an  activity  drop  below  the  10  “  activity  threshold  at 

approximately  31.8  years  ( 109  s)  after  fuel  irradiation.  For  times  longer  than  10.0 
seconds,  the  curves  from  each  input  problem  closely  resemble  each  other  with  the 

exception  of  some  separation  around  the  one  year  (3.15*10  s)  point.  This  separation  is 
more  clearly  shown  with  a  semi-log  plot  of  the  data  as  in  Figure  22.  The  straight,  parallel 
lines  are  a  reflection  that  all  of  the  activity  at  long  times  is  caused  by  the  same  isotope, 
Rh-106  ( thalj-  =  28.9s).  Rh-106  is  in  secular  equilibrium  with  its  parent  isotope,  Ru-106 

( half  =1-02  years).  There  is  a  99.9%  correlation  between  the  total  quantity  of  initial 

fission  products  with  a  mass  number  equal  to  106  and  the  ratio  of  the  activities  shown  in 
Figure  22. 
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Figure  20  Gamma  Activity,  3-4  MeV,  1  ns  to  10  s,  Log-Linear 
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Figure  21  Gamma  Activity,  3-4  MeV,  1  ns  to  100  y,  Log-Log 
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Figure  22  Gamma  Activity,  3-4  MeV,  1  ns  to  100  y,  Semi-Log 
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5.  4-5  MeV 


The  activity  due  to  gamma  radiations  with  energies  between  4  MeV  and  5  MeV 
continues  the  trend  of  decreasing  contribution  to  the  total  activity,  as  is  shown  in  Figure 
23.  The  crossing  of  some  of  the  activity  curves  is  shown  clearly  in  Figure  24  in  a  semi¬ 
log  plot.  Crossing  activity  curves  for  a  fission  product  decay  problem  indicates  that  the 
new  activity  comes  from  not  just  a  new  set  of  isotopes,  but  that  the  activity  from  a  new 
mass  chain  has  become  predominant. 

Figure  25  shows  the  late  time  activity  of  the  different  fuel  sources.  At 

approximately  5.3  days  ( 4.6  *  105  s)  after  irradiation  with  high  energy  neutrons,  the  Pu- 
239  source  shows  a  departure  from  the  curves  of  the  other  seven  test  scenarios. 

This  effect  is  shown  more  prominently  on  a  semi-log  plot  as  in  Figure  26.  The 
different  slopes  indicate  a  difference  in  isotopes  that  contribute  to  the  activities. 
Inspection  of  the  data  show  that  the  activity  curves  of  all  of  the  isotopes,  except  Pu-239 
high,  are  more  than  99.99%  made  by  the  decays  of  Rb-88  ( thalj-  =  17.73  minutes).  The 

relatively  short  half  life  of  Rb-88  accounts  for  the  very  low  activities  present. 

Fission  products  are  typically  neutron  heavy  isotopes  that  under  go  beta  decays. 
However,  the  fission  of  a  Pu-239  nucleus  with  a  high  energy  neutron  results  in  the 
production  of  Ga-66  ( thalj-  =  9.49  hours),  a  neutron  light  isotope  that  undergoes  an 

electron  capture  decay.  Ga-66  is  only  produced  on  average  at  a  rate  of  2.02  *  10~14  atoms 
per  fission  of  Pu-239  by  a  high  energy  neutron  but  produces  enough  high  energy  gammas 
to  be  significant  as  the  activity  from  other  fission  products  decays  away. 
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Figure  23  Gamma  Activity,  4-5  MeV,  1  ns  to  10  s,  Log-Linear 
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Figure  24  Gamma  Activity,  4-5  MeV,  1  ns  to  10  s,  Semi-Log 
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Figure  25  Gamma  Activity,  4-5  MeV,  10  s  to  106  s,  Log-Log 
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Figure  26  Gamma  Activity,  4-5  MeV,  1  s  to  1.2*106  s,  Semi-Log 
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6.  5-6  MeV 


The  activities  produced  from  fission  in  the  5-6  MeV  range  differ  significantly 
from  previously  discussed  energy  bins.  The  different  test  cases  are  separated  by  the  fuel 
type  rather  than  incident  neutron  energy  but  U-235  now  exhibits  a  greater  activity  than  U- 
238.  Pu-239  activity  is  still  the  lowest  of  any  fuel.  Figure  27  also  shows  an  increase  the 
activity  for  the  U-238  fuels  between  1.0  second  and  10.0  seconds.  Figure  28  shows  this 
time  scale  in  semi-log  fonnat  to  highlight  the  differences  in  the  curves.  The  top  five 
gamma  intensity  contributions  to  the  activity  at  1 .0  second  are  all  the  result  of  the  decay 
of  Rb-92  ( thalj-  =  4.49s).  The  same  analysis  at  3.0  seconds  yields  the  same  result;  the 

decay  of  Rb-92  is  still  responsible  for  the  top  five  contributors  to  the  total  activity  in  the 
5-6  MeV  bin. 
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Figure  27  Gamma  Activity,  5-6  MeV,  1  ns  to  10  s,  Log-Linear 
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Figure  28  Gamma  Activity,  5-6  MeV,  1  ns  to  10  s,  Semi-Log 
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The  cause  for  the  spike  in  activity  in  U-238  fast,  and  to  a  lesser  extent  U-238 
high,  is  actually  a  result  of  the  initial  distribution  of  fission  products.  While  relatively 


little  Rb-92  is  present  initially  in  U-238  compared  to  U-235,  the  precursors  to  this  isotope 
(only  the  direct  decay  branches  are  shown) 


92  Sp 

34ae- 


pr 


3.00*10"',? 


+  35Br 


P  ,  92^„  92R, 

0.343?  >  36  1.840?  >  37  b 


are  present  in  greater  quantities  in  U-238.  The  initial  quantities  as  shown  in 
Table  5  and  the  results  after  3.0  seconds  are  shown  in 


Table  6. 


Table  5  Fission  Product  Yield  per  Fission  Type 


Isotope 

Neutron 

Energy 

Se-92  (t=0) 

Br-92  (t=0) 

Kr-92  (t=0) 

U-235 

thermal 

4.17E-07 

2.68E-04 

1.66E-02 

fast 

5.51E-07 

2.09E-04 

1.35E-02 

high 

1.15E-06 

3.03E-04 

8.25E-03 

U-238 

fast 

1.98E-05 

2.34E-03 

2.50E-02 

high 

9.91E-06 

1.34E-03 

1.58E-02 

Table  6  Rb-92  Quantity  per  Fission  at  t=0  s  and  t=3  s 


Isotope 

Neutron 

Energy 

Rb-92 

(t=0) 

Rb-92 

(t=3) 

U-235 

thermal 

3.13E-02 

2.85E-02 

fast 

2.79E-02 

2.49E-02 

high 

2.93E-02 

2.29E-02 

U-238 

fast 

1.47E-02 

2.33E-02 

high 

1.89E-02 

2.07E-02 
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Figure  29  shows  an  even  more  defined  separation  of  Pu-239  activity  following 
fission  by  a  high  energy  neutron  than  demonstrated  in  the  4-5  MeV  bin.  The  explanation 
for  the  difference  is  the  same  as  above,  namely  the  presence  of  Ga-66.  Here,  the  only 
difference  is  that  the  activity  in  the  other  test  cases  is  caused  by  the  decay  of  Rb-90  as 
opposed  to  Rb-88  in  the  4-5  MeV  case. 
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Figure  29  Gamma  Activity,  5-6  MeV,  10  s  to  106  s,  Log-Log 
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7.  6-7  MeV 


The  6-7  MeV  energy  bin  activity  displays  many  similarities  to  the  5-6  MeV  bin. 
In  fact,  the  activity  at  early  times  is  due  to  the  same  isotope,  Rb-92.  The  rest  of  the 
curve,  as  shown  in  Figure  31,  depicts  the  activities  from  each  of  the  different  sources  as 
largely  indistinguishable.  The  semi-log  plot  in  Figure  32  is  more  interesting  because  of 

—i  ? 

the  t  '  plot  as  well  as  the  closeness  and  straightness  of  the  activity  curves.  The  semi¬ 
log  plot  clearly  shows  that  in  the  6-7  MeV  energy  bin,  the  approximation  is  a  poor  fit  for 
the  data.  The  tight  packing  of  the  curves  comes  from  all  of  the  top  live  activity 
contributors  being  products  of  1-136  ( thalf  =  83.4  s)  decay. 
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Figure  30  Gamma  Activity,  6-7  MeV,  1  ns  to  10  s,  Log-Linear 
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Figure  31  Gamma  Activity,  6-7  MeV,  10  s  to  10,000  s,  Log-Log 
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Figure  32  Gamma  Activity,  6-7  MeV,  10  s  to  6,000  s,  Semi-Log 
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8.  7-8  MeV 


The  activity  in  the  7-8  MeV  bin  is  unique  in  that  for  the  test  problems  studied,  the 
activity  in  this  bin  is  the  result  of  only  one  gamma  ray  from  one  isotope,  the  7.000  MeV 
gamma  emitted  by  the  beta  decay  of  Br-88  (thay  =  16.29  s). 

The  increase  in  activity  after  1.0  second  is  due  to  the  fission  product  yield.  The 
fission  product  decays  that  result  in  the  formation  of  Br-88  are  shown  below  (only  the 
direct  decay  chain  is  shown) 


Se-89  ( thaij-  =  0.41  s)  also  results  in  some  production  of  Br-89  through  a  beta  decay 


followed  by  a  neutron  emission.  The  combination  of  the  A=88  chain  with  the  Se-89 
decay  accounts  for  the  increase  in  the  activity  of  the  U-235  and  U-238  curves. 

Figure  34  demonstrates  again  the  application  of  the  t  '  approximation. 
The  semi-log  plot  shows  that  the  decay  continues  to  be  that  of  a  single  isotope,  Br-88. 
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Figure  33  Gamma  Activity,  7-8  MeV,  1  ns  to  10  s,  Log-Linear 
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Figure  34  Gamma  Activity,  7-8  MeV,  1  s  to  1,200  s,  Semi-Log 
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9.  >8  MeV 


Energy  bins  with  minimum  energies  above  8.0  MeV  do  not  have  any 
activity  for  the  test  problems  examined.  Gamma  radiation  energies  greater  than  this 
threshold  are  due  to  the  decays  of  low  Z  materials  (C-15,  Na-20,  Al-24,  P-28,  and  K-36), 
which  have  a  negligible  probability  of  being  a  fission  product. 
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VII.  Conclusions  and  Recommendations 


VII.  A :  Conclusions 

The  primary  objective  of  this  thesis  was  to  produce  a  more  precise  code  for 
calculating  complex  decay  chains  with  greater  precision  than  available  by  current 
methods.  In  implementing  the  exponential  moments  function  with  the  use  of  a 
transmutation  matrix,  a  code  has  been  developed  which  gives  the  user  access  to  quick  and 
precise  results  for  problems  of  as  much  complexity  as  desired.  The  code  also  provides 
insight  into  the  effect  on  activity  caused  by  the  uncertainty  of  isotope  data. 

In  implementing  the  exponential  moments  function,  a  proven  mathematical  model 
used  in  neutron  transport  calculations,  the  data  was  manipulated  to  provide  the  correct 
arguments  necessary  to  calculate  decay  values.  Verification  of  the  method  was 
performed  by  comparison  to  a  Mathematica  benchmark  implementing  the  Bateman 
solution  formula.  The  strong  correlation  between  the  two  calculation  methods  gives 
strong  indication  that  the  exponential  moments  function  is  in  fact  accurate  out  to  10  digits 
when  using  double  precision  values  in  the  calculation. 

The  use  of  a  transmutation  matrix  is  a  novel  way  of  approaching  the  decay  chain 
problem.  The  matrices  are  stored  in  a  sparse  format  which  means  that  a  library  can  be 
constructed  to  significantly  reduce  calculation  time  for  some  problem  sets.  The  primary 
advantage  of  the  T-matrix  is  that  it  allows  the  code  to  evaluate  many  different  decay 
problems  simultaneously. 

This  thesis  also  demonstrated  one  possible  application  of  the  exponential 
moments  function  combined  with  the  use  of  a  T-matrix.  This  approach  proved  a 
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worthwhile  method  that  lead  to  an  increased  understanding  of  the  activity  emitted  by 
fission  products  with  respect  to  both  energy  and  time. 

The  overall  results  of  the  thesis  is  that  a  method  for  calculating  decay  products 
with  increased  precision  has  been  developed  and  can  be  adapted  by  the  community  where 
such  a  requirement  is  needed  to  better  understand  complex  decay  problems  or  where  the 
uncertainty  of  gamma  activity  or  isotope  quantity  is  of  interest. 

VII.  B:  Recommendations  for  Future  Work 

To  implement  the  exponential  moments  function  calculations  and  generate  T- 
matrices,  certain  assumptions  had  to  be  made  to  ensure  that  a  working  code  could  be 
made  and  thoroughly  verified.  Further  exploration  of  some  of  these  assumptions  could 
produce  even  better  results.  Daughters  of  spontaneous  fission  decays  were  neglected  in 
this  code.  Also,  the  treatments  of  metastable  states  in  this  thesis  were  not  complete. 
Future  work  aimed  at  adding  spontaneous  fission  decays  and  adding  more  realistic 
metastable  state  tracking  and  calculation  should  have  the  largest  impact  in  improving 
performance  even  farther. 

Further  work  in  optimizing  the  exponential  moments  function  to  run  faster  would 
also  be  a  benefit  to  the  code.  In  addition,  this  code  used  the  built  in  pseudo-random 
number  generator  that  is  part  of  Fortran  90/95.  The  implementation  of  a  better  pseudo¬ 
random  number  generator  and  performing  more  than  1000  runs  would  improve  the  results 
of  the  Monte  Carlo  simulation. 
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Appendix  A:  Order  Invariance  of  Bateman  Equation 


The  Bateman  equation  for  the  kth  isotope  is  given  as 


f  k—\  \  k  —  Xjt 

K(o=N,(  o)  nw  — 

V'=1  J 


(A.l) 


I  start  the  analysis  by  separating  out  two  arbitrary  isotopes  from  the  decay  chain,  m  and  n. 


Nk(t)  =  Nt(0)bmm+]bnn+]AmtAnt 


k- 1 


nw.' 


i= 1 
i^m 
y  i*n 


-A,t 


-  e  ' _ + _ LI _ + _ e— _ 

(An-Am)mAlt-Amt)  ( Am-An)U(Ait-A„t ) 


Jrm  i= i 
J*n  i*j 


i=i 

i^m 


i= 1 
i^n 


(A.  2) 


The  next  step  is  to  swap  m  and  n. 


Nk  (t)  =  Aj  (0)b  +1b  +lAntAmt 


f  \ 

rk^ 

i=i 

i^n 


-At 

e  " 


Zt - + - 

(z„-z,) (A„-A„)n (Z,t-Zj) 


l*J 


(A. 3) 


Because  equations  (A.2)  and  (A. 3)  are  equal,  the  Bateman  equation  is  invariant  with 
respect  to  order  of  the  argument. 
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Appendix  B:  Depth-First  Search  Verification  Test  Problems 

1.  A^B 

„  A^B 

2. 

-^C 

A^B^-E 

-^C^D^E 

A^B 

4.  ->C 
— »  D 

A^B^D 

5.  -^B^E 
-^C 

,  A^B^D 

6. 

-^C^D 

A^B^D 
->  C  ->  E 

8.  A^B^C^D^E 

9.  A^B^C^---^T^U 
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Appendix  C:  Mathematica  Verification  Notebook 


Off [General: :spelll] ; 

SetDirectory  [  "C :  \Thesis\ Verify"  ]  ; 
prdblemData  =  Import [  "TestLanfcrias . csv" ,  "CSV"]  ; 
imax  =  Length  [problenData] 
strm  =  CpenWrite  [  "answers .  txt"  ]  ; 


starttime  =  AbsoluteTime  [  ]  ; 

Do[ 

prcblarData  [  [  i ,  1]  ] 

A[l]  =  N[10  KOI  ,  300]  ; 

prcblsrData  [  [  i ,  2]  ] 

A [2]  =  N[10  K®  ,  300]  ; 

paxblaTData[  [i,3]  ] 

A [3]  =  N[10  1003  ,  300]  ; 

A[4]  =  N[l0P^l«Data[[i,4]]  ,  ^  . 

D°[ 


bateman[k_] 


-A[j] 


(ni=i  (^[i]  -  Atj]))  oifcj+i  <*[*■]  -  Atj])) 


b  =  bateman[k]  ; 

Write[strm,  FortranForm[N[b,  20]  ]  ]  ; 
,  {k,  4}] 

t  {ij  imax}] 

steptime  =  AbsoluteTime  [  ]  ; 

Close  [strm]  ; 

elapsedtime  =  stoptime  -  starttime 
Quit[]  ; 
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